Report on Thermal Neutron Diffusion Length 
Measurement in Reactor Grade Graphite Using 
MCNP and COMSOL Multiphysics 



School of Physics and Astronomy University of Birmingham Birmingham B15 2TT United 
Kingdom 

E-mail: srml05@bhain.ac.uk 

Abstract. Neutron diffusion length in reactor grade graphite is measured both experimentally 
and theoretically. The experimental work includes Monte Carlo (MC) coding using 'MCNP' 
and Finite Element Analysis (FEA) coding suing 'COMSOL Multiphysics' and Matlab. 
The MCNP code is adopted to simulate the thermal neutron diffusion length in a reactor 
moderator of 2m x 2m with slightly enriched uranium (^^^U), accompanied with a model 
designed for thermal hydraulic analysis using point kinetic equations, based on partial and 
ordinary differential equation. The theoretical work includes numerical approximation methods 
including transcendental technique to illustrate the iteration process with the FEA method. 
Finally collision density of thermal neutron in graphite is measured, also specific heat relation 
dependability of collision density is also calculated theoretically, the thermal neutron diffusion 
length in graphite is evaluated at 50.85 ± 0.3cm using COMSOL Multiphysics and 50.95 ± 0.5cm 
using MCNP. Finally the total neutron cross-section is derived using FEA in an inverse iteration 
form. 



1. Introduction 

This work demonstrates an analytic approach accompanied with models of Finite Element 
Analysis (FEA) and Monte Carlo (MC) with an experimental measure on neutron cross-section 
and slowing down process. In MC approach Monte Carlo N-Particle Transport Code (MCNP) 
is used to simulate the simplified version of reactor moderation process. Similarly in FEA the 
moderator modelled (Assuming a symmetrical distribution) using point kinetic equations, based 
on partial and ordinary differential equation in software package. 

2. Theoretical Calculations 

Having the number of particles found in a volume element dr where dr = dxdydz at r with a 
vector with solid angle at Q. be donated by [Ij: 
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N{r,n,t)drdn 




Therefore can have: 



dN 
~dt 

dN 



-NVa + J N{r, Q' ,t)Vasf{n.n')dQ' + S{r,n,t) 



(2) 



Where the first term donates the number of particles present in given volume (particle den- 
sity) and second term —NVa represents the total number of particles removed from the given 
volume by scattering and capture, a is representing the total cross-section. The third term 
represents the total number of particles scattered into the given volume, and f{Q.il') represents 
the relative probability of scattering through an angle whose cosine is Jl.Q', where is a unit 
vector in the direction of the initial velocity and is unit vector in the final direction. Finally 
S{r, Q, t) is the external source term available in the system and is given by: 



N{r, n, t)drdn = N{Z, if)dZdipd<j) 



(3) 



Where (p is the cosine of the velocity vector in the Z direction and </) is the longitude of velocity 
vector. 




Figure 1. The Velocity Vector: Where ip is the cosine of the velocity vector in the Z direction 
and (p is the longitude of velocity vector and is unit vector in the final direction. 



Now an assumption can be made such: 
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Rewrite the equation [T] as: 



5N{Z,p) 



No{Z) = 2iT I N{Z,ip)dip 



-N{Z, ip)Va + J N{Z, ip')Vasf{ipo)dn + Siz) 



(4) 



(5) 



^ dZ 

Where (fo is the cosine of the angle between initial and final velocities and it can be found by: 



cosOcosO' + sin9sin9'cos{(j) — cp') 



(6) 



Or using: 

ipo = ipif' + ^Jl-Lp^^Jl- ip'^cos{4> - (j)') (7) 

Now if the collision function of F(lpq) expanded in spherical harmonics: 

°° 2/ -I- 1 

^('^o) = E^^i^i('^o) (8) 



With Fi = J f{ipQ)Pi{ifQ)dQ. Using similar expansion the phase density function will be: 

NiZ,^) = Y.^N,{Z)P^{^) (9) 



Where N{Z,(p) = J N[Z,ip)Pi[ip)dQ,, with assumption that N{Z,ip) is isotropic, three condi- 
tions must be satisfied all the times: first, it is far from the source (equal to Mean Free Path 
(MFP)), second, it is far from the boundaries; third, the probability of capture is small compared 
to probability of scattering. Having all the conditions satisfied, the following can be assumed: 

N{Z,^)^^{Nq{Z) + 3^Ni{Z)) (10) 

Where the second term in the bracket donates the particle flux (J). Here Ni = J LpN{Z, ip)d^l = 
J/V. For simplicity we choose our unit such that V = 1 and a = 1, hence: 

1-/ = - (11) 
a 

Now the Boltzmann equation takes the form of: 

dN 



N + {l-f)l N(Z,v')/(No)dn' + S{Z) (12) 



By integrating the equation over all possible angles (dO) we have: 

dN^ 



dZ 



i = -A^o + (1 - f)No + 4ttS{Z) (13) 



/ is normalized in such a way that / f{^}.n')dn' = J f{n.n)dn = l, hence going back to Eq. [o] 
for the case / = we have: 

Fo = J f{^o)dn = 1 (14) 

Hence by integrating over all angles and Multiplying by ip we have: 

l^ = -N^ + {l-mN^ (15) 
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Now the second order differential equation gives: 



3(1 - /i dZ^ 
This also can be written as: 



^ -fNo + S,{Z) (16) 



iVo - ^A^o + ^^0 = (17) 



Eqs. 16 and 17 are known as diffusion equation. Here L is diffusion Length abd D is diffusion 
coefficient and it is equal to gjzV = , where and Xtr are the scattering and transport 
mean free path. Xtr can be calculated from: 

^tr = . (18) 
1 - {C0s9)av 

and {cos6)av is equal to / f {ipQ)LpQdVL = /i. Also can be measured using following relation: 

= ^ (19) 

Here Ac is the capture mean free path. 
MAXIMUM ENERGY LOSS 

If a neutron with initial velocity Vq collides with a nucleus of mass M (at rest), then in the 
Centre of Mass (CoM) system, the initial velocity is MVq/ M + \ after collision. The momentum 
of of neutron and the nucleus will be equal to magnitude oppositely directed vector. Figure [2] 
demonstrates the collision in CoM system. 

As demonstrated in Fig. [2] the 9 is the deflection angle and G is angle on the final velocity v. 
The f in this case is given by: 

cosO H = vcos'a (20) 

M + 1 M + 1 ^ ^ 



so 



since u = log%- then: 



+ - ^cose = (21) 

w + r w + r M + 1 ^ ' 



cosO = l-^^^l-^f (22) 



(M + 1)2 

cos9 = 1 - ^ ^ ' 1 - e-" (23) 
2M ^ ' 
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Figure 2. The Collision in Centre-of-Mass System:If a neutron with initial velocity Vq collides 
with a nucleus of mass M at rest, then in the Centre of Mass (CoM) system the initial velocity 
is MVq /M + 1 after collision. The momentum of of neutron and the nucleus will be equal to 
magnitude oppositely directed vector. Here 6 is the deflection angle and B is angle on the final 
velocity v. 



now differential cross-section gives: 

dcos9 _ _ (M + 1)2 
du ~ 2M ' 

Hence: 



(24) 



(M + l)2 ^ M-1 ^ , ^, 

COS@ = —- ^6 2 6 2 (25) 

Therefore the maximum logarithmic energy loss can be calculated from: 

QM = iogij^f (26) 

The qm is at most when Q = n. Now going back to the problem we can redefine the collision 
density function as: 

The term ^^^^ e~" is the normalization constant chosen to satisfy / dO. J duf{ipo,u) = 1 and 
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6 is the Dirac 6 function. So that 6{x — a) = when x ^ and / d{x — a)F[x)dx = F(a). Now 
the average logarithmic loss (^) can be calculated from: 



e = 1 - 4^^".^-- (28) 



and {cos)aV is 2/3M. 



Energy Distribution of Slowed Down Neutrons 
I. Stationary Case 

The average collision density per unit time, with logarithmic energy intervals is given by: 

^q{u) = r du'^Q{u')h{u')h{u - u') + 5{u) (29) 







where /o(ii) is (M + 1)^6 "/4m for u < qm and it is zero otherwise. In stationary case the total 
number of neutron produced is unity per unit in this case, i.e. for M = 12, / + 0{u) becomes 



3.5e Hence the equation 29 becomes: 



Mu)= dMVo(n')/i(^t')3-5e-("-"') + 5(m) (30) 

J u—qm 



where qm = 0.72 



II. Time-dependent Case 

The time dependent when there is no absorption in the system and source strength is unity and 
is given by: 



l{u) dipQ 

~~6r 



+ iPo{u,t)= du'Mu,t)e-'^''-'''^ +6{u)6{t) (31) 







where l{u) is the mean free path and if the mean free path is constant, the Laplacian form of 
the equation for M ^ 1 becomes: 

s/n r'^ 

l + ^Mu,s)= du'Mu',s)fo{u-u') + 6{u) (32) 



V Jo 



now: 



w 



A( \ ~ I ^ ,w'(p{w',s) 



Eq. 33 applies for u > qm where w = Iqs/v, r = M — 1/M + l and 4>{w, s) = (1 + w)4>o{u, s), so 



that the mean free path is proportional to velocity. 
III. Rigorous Numerical Solution 

The slowing down process is not an easy approach, therefore a more discrete form of solution 
also could be defined using: 

F{E) = / ^(^' ^ E)(l){E')dE' + 5{u) (34) 



Where J2si'^' ~^ ^) scattering term between energies E' and E. Recalling 

Y^{E' ^E) = Y^{E')P{E' ^E) (35) 

s s 

Where P{E' — )■ E) is the probability of collision happens between E' and E. It can be defined by: 

P{E' ^ E).E'{1- a) = 1 (36) 



Hence: 



Now: 



For M = 1 the collision is defined as: 

Es(-B') 



^ ^' = (1^ 



F(E) = 1^ ^^,KE)dE' + 4(E) (39) 



Also the solution with capture process: 



Also For M 7^ 1 the collision density can be found: 

F{E) = ^dn^,^{E')dE' + (41) 
JE/a (1 - a)E' (1 - a)Eo 

This only applicable if aE^ < E < Eq. A theoretical calculation is performed for an arbitrary 
system and Fig [3] is derived. For graphite the collision density is also measured for different 
neutron energy range as demonstrated by Fig. |4| 

and for graphite: 

The oscillations are due to Plaezack Oscillations which is the fundamental phenomenon 
associated with the neutron slowing-down [2j. And finally for the M ^ 1 with capture: 



F{E) = {J2{E) + J2{EmE) 



^ EsiE'mE') ^^, So \ 

E/a {l-a)E' {l-a)EoJ 

So 



JE/a EtiE') (1 - a)E' ^ (1 - a)Eo 
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3. Model Set-up in MCNP 

The MCNP code is developed in Los Alamos National Laboratory and it is well-known for 
analysing the transport of neutron and 7-rays in matter. MCNP is a continuous energy mod- 
eller with generalized geometry time dependent code that implements data from nuclear libraries 
such as, Evaluated Nuclear Data File (ENDF), Evaluated Nuclear Data Library (ENDL), Acti- 
vation Library (ACTL). 

The code structure is divided into four main sections, geometry definitions, surface definitions, 
material cards, and tallies. Geometry of MCNP is a three-dimensional form defined using cell 
and surface cards. For instance Fig [5] demonstrates the geometry setup in this system. FigurejH]- 
illustrates the reactor moderator, where each cylinder represents the fuel rod containing slightly 
enriched ^^^[/. The moderator is reactor grade graphite. 

The user can instruct the code to make various analysis with tally cards. The tallies are to 
measure the particle current on the surface, particle flux and energy deposition. In fact any 
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Figure 5. The MCNP Geometry Set up: Figure demonstrates the reactor moderator module. 
The dimension was set 2m x 2m 



quantity in form of Eq. 43 can be talhed [3] . 

C = j ^{E)f{E)dE (43) 

Here (j){E) represent the particle flux and f{E) is the cross-section quantities given in the li- 
braries. Table [T] demonstrates the six MCNP standard tallies. 
In MCNP when neutron collides with a nucleus: the nuclide will be identified depends on the 



Table 1. MCNP Neutron Tallies 

Property Data 

F1:N Surface Current 

F2:N Surface Flux 

F4:N Track Length Estimate of Cell Flux 

F5a:N Flux at a Point 

F6:N Track Length Estimate of Energy Dependence 

F7:N Track Length Estimate of Fission Energy Dependence 



preferences of target, that is either the S{a^[i) treatment or velocity of target; therefore the 
nucleus will be sampled for low energy neutrons: neutron capture or absorption will be modelled 



and either elastic or inelastic reaction depend on the model performance. 

However sometimes different nuclide form a material, (where the collision occurs) therefore we 
can have: 

k—l n k 

EE<eEE<EE m 

1=1 ti 1=1 ti 1=1 ti 

Where X)tj is the microscopic total cross-section of nuclide i. The total cross-section is sum of 
the capture cross-sections in the cross-section reference table. 

The collision between thermal neutrons and the target will be effected by thermal motion of the 
atoms, chemical binding and lattice structure of the target. This is called Free Gas Thermal 
Treatment. Hence the effective scattering cross-section in laboratory system is given by [4J: 

alf^{E) = -11 as{vref)vreiP{v)dv^ (45) 

Here Vn is particle velocity, Vrei is the relative velocity, P{v) is the probability density function 
and 9? as explained before is the cosine angle of velocity vector. The relative velocity can there- 
fore is given by: 

Vrel = {v^n + V^ -^VnV^t)^^^ (46) 



The density function is also given by: 



P{v) = -^l3-^vh-^\^ (47) 



where f3 = (41^)^^^- However most of the time in equation 45 the ag can be ignored for heavy 
nuclei, where arei can have moderating effect and is given by|4j: 



P{v, if) oc \lvlv'^ - 2?;.t;„(^it;2e-/3'^' (48) 

In MCNP there are also two types of capture, analogue and implicit. Analogue occurs when 
the particle is killed with probability of Where aa and at is the absorption and total cross- 
section respectively. Implicit capture happens when neutron weight iWn) is reduced by number 
of collisions and is given by: 

Wn+l = (1 - —)Wn (49) 

The elastic scattering directed by two body kinematics: 

Eout = \Ein{{l - a)ecm + 1 + a) (50) 

Where is the center of mass cosine of angle between incident and existing path direction. 
Where in inelastic an scatter the particle reaction is chosen such as (n, n'), (n, 2n), (n,/), and 
(n, n'a), and is given by jl|: 

,,^,.^^^ E^2B(A.r,^ ^ (51) 
10 + ^1 



and 




Vlab = @cm\^ + ^r—\ -FJ (52) 



Here ipiab is cosine of laboratory scattering angle. However for thermal energy neutron, 
S{a, (3) treatment is needed. For inelastic treatment the secondary particle distribution will 
be represented by set of discrete energies between 4ey to \{)~^eV . 

4. Model Set-up in COMSOL 

The Partial Differential Equation (PDE) module of COMSOL package supports three types of 
formation: coefficient form, general form, weak form. The coefficient form is a linear system 
where as the general and weak form supports non-linear, and more flexible form of definitions 
is supported by weak form. In this report one study is performed for thermal group transport 
using equation based general form of the system. 

In equation based system the independent variable ui, U2 will be defined in following equation: 
S'^u Su 

'^'^"■'St ~ '^'^'^ V u + au - j) + (3. \j u + au = f (53) 

Where is the mass matrix, and eajjj is called mass term, daff is called damping term, 
\J .{c u + au — ^) is called diffusive flux, f3 is convection flux, a is the absorption coefficient, 
and / is the source term. 

Environmental factors are defined by enforcing boundary conditions using Dirichlet equation. 
Dirichlet imposes Laplace equation (our transport equation) to the system domain. It is there- 
fore more convenient to have the numerical Laplace such that: 

5^U 5U 

Now that Laplace equation is defined we need to numerically define the flux and multiply the 
two values, hence: 

\/ .{-c\/ u- au + 'y) (55) 



Equation [55] is called flux vector. Here a is the velocity term, 7 is the source term, c can be 
also indirectly calculated for an anisotropic material. 

Equation 53 is in computational domain (a;), thus the calculations need to satisfy all conditions 



in boundary domain. This is called Neumann-Dirichlet where the boundary will be transformed 
from uj to duj (from computational boundary to domain boundary). This transformation is also 
described as domain decomposition preconditioner [5]. Thus the partial differential equation is 
given by + u = where y donated as Laplacian therefore we will have: 



du 
dn 



x) = \/'^u{x).n{x) 



(56) 



Where n refers to a normal vector, thus we can rewrite the equation 55 as: 



n.{c u + au — ^) = g — /i^ 
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(57) 



where g and h^j^ donate the boundary source term and the Lagrange multipher factor, /i^ is 
needed in a mixed field situation as it corresponds to local maxima and minima. In some respect 
/i^Y can also refers to the velocity [6j. 

Now by taking the energy dependent diffusion equations we have: 



-^.Vi^V0 + E'^= rj2(E^E')<p{E')dE'+x{E) [ v{E')Y.iE')HE)d{E') + S{r,E,t) 

V Ot f Jo s f 

(58) 

Where in multi-group theory discrete energies varies with G discrete group as: 
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Figure 6. Discrete Group Relation, Thermal, Epithermal, Fast Region 



The group flux can be obtained by integrating total fluxes across the group energy range. Hence 
the parameters can be defined as below: 
I. Total Cross Section: 



g 

t 



Y^{E)<P{E)dE 



j:tiE)<i>{E) 



dE 



II. Diffusion Length: 



III. Inverse Velocity: 



Jj,.EgD{E)vct>{E)dE 



lE,-iE9VHE)dE 

1 _ lE,-lEg^)Cj^{E)dE 



(59) 



(60) 



(61) 



IV. Fissile Spectrum Term 



Xg 



Egx{E)dE 



(62) 



Therefore the stationary solution for many group equation can be given by(in this work the 
equation is only solved for thermal group spectrum): 



1^ 

Vg 5t 



G g'^g g' 
* g'=\2 " 9-1 / 



(63) 



Where the right and left hand side of the equation present the loss and production term 
respectively, v is the average neutron speed Xg is the fraction of prompt neutrons. For the 
simulation the Arbitrary Lagrangian Eulerian (ALE) mapping mesh analysis is used. 

5. Discussion and Results 

This report has reviewed the neutron diffusion length both using COMSOL Multiphysics and 
MCNP. The total number of 5,000,000 meshes used for the iteration process in COMSOL. 
Both the thermal neutron flux and absorption property of graphite with respect to its cross- 
section features have been evaluated. The thermal diffusion length therefore calculated was 
50.85 lb 0.3cm in COMSOL and 50.95 it 0.5cm in MCNP. Figure [7| demonstrates the distribution 
of thermal neutron increases as they penetrate deeper into graphite compared both in MCNP 
and COMSOL. The red line in the figure is also illustrates the experiment done in the lab on 
graphite using Am — Be source. The Am — Be was canned on top of an aluminium cylindrical 
tube. Two set-up is used in this experiment, a cadmium cover with nominal thickness of 0.1 
cm (As explained previously the cadmium has cut-off of 0.55ey) and were constructed to fully 
overlap the detector edges to avoid leakages. The flux distribution is measured by putting the 
source at a fixed location and relocating the detector at 25 cm distance intervals in horizontal 
and vertical directions. 
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Figure 7. Figure demonstrates the the neutron diffusion length using COMSOL and MCNP, 
the circles illustrates the COMSOL results where as the stars demonstrate MCNP. The red line 
as well shows the experiment was done using BF3 tube. 
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As shown in figure [7] the distribution calculated using COMSOL is less than half order of 
magnitude higher than MCNP. For completion the absorption probability cross-section in FEA 
is also evaluated using inverse iteration technique as demonstrated by figure |9j It is clearly 
shown as the neutron travels deeper into graphite they probability of absorption in graphite is 
also increases. 

To understand the respond of absorption cross-section to different thermal neutron energies, 




Figure 8. The Total Neutron Absorption in Graphite: Derived Using COMSOL. It shows the 
probability of absorption increases as the neutrons travel deeper in to the graphite, the areas 
red shows the highest and blue the lowest. 

the evaluated values are compared with with j7] , [8] , [9] , [lO] , [H] , [l2] and [l3j as demonstrated 
in figure |9| 
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Figure 9. The Total Neutron Absorption Cross-section in Graphite Compared with [7], [S], [S]; 
ilOj) iUj) and Fittings are extracted by permission from [14] 
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